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Abstract 

We  have  previously  developed  a phase-field  model  of  solidification  that  includes  con- 
vection in  the  melt  [Anderson,  McFadden  and  Wheeler,  Physica  D,  135  (2000)  175-194]. 
This  model  represents  the  two  phases  as  viscous  liquids,  where  the  putative  solid  phase 
has  a viscosity  much  larger  than  the  liquid  phase.  The  object  of  this  paper  is  to  examine 
in  detail  a simplified  version  of  the  governing  equations  for  this  phase-field  model  in  the 
sharp- interface  limit  and  derive  the  interfacial  conditions  of  the  associated  free-boundary 
problem.  The  importance  of  this  analysis  is  that  it  reveals  the  underlying  physical  mech- 
anisms built  into  the  phase-field  model  in  the  context  of  a free-boundaxy  problem  and,  in 
turn,  provides  a further  validation  of  the  model.  In  equilibrium  we  recover  the  standard 
interfacial  conditions  including  the  Young-Laplace  and  Clausius-Clapeyron  equations  that 
relate  the  temperature  to  the  pressures  in  the  two  bulk  phases,  the  interface  curvature  and 
material  parameters.  In  nonequilibrium  we  identify  boundary  conditions  associated  with 
classical  hydrodynamics,  such  as  the  normal  mass  flux  condition,  the  no-slip  condition  and 
stress  balances.  We  also  identify  the  heat  flux  balance  condition  which  is  modified  to  ac- 
count for  the  flow,  interface  curvature  and  density  difference  between  the  bulk  phases.  The 
interface  temperature  satisfies  a nonequilibrium  version  of  the  Clausius-Clapeyron  relation 
which  includes  the  effects  of  curvature,  attachment  kinetics  and  viscous  dissipation. 


PACS:  05.70.Fh;  05.70.Ln;  81.10.Fq 
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1.  Introduction 


In  a recent  paper  [1]  (AMW)  we  derived  a phase-field  model  that  describes  the  solidification  of 
a pure  material  when  fluid  motion  is  present.  The  motivation  for  developing  and  studying  such 
a model  follows  from  experimental  observations  and  theoretical  analyses  which  show  that  fluid 
motion  in  the  melt  can  have  a profound  effect  on  the  morphology  of  the  solid-liquid  interface  [2]— 
[6].  In  an  industrial  setting  morphological  and  hydrodynamic  instabilities  can  have  a detrimental 
effect  on  both  the  control  of  the  process  and  the  quality  of  the  solidified  product.  However,  even 
when  careful  attempts  are  made  to  suppress  fluid  motion  by  eliminating  external  forces  such 
as  gravity,  it  is  not  always  possible  to  eliminate  all  fluid  motion.  For  example,  expansion  or 
contraction  flow  is  present  during  solidification  when  the  solid  and  liquid  densities  of  the  material 
are  different.  Furthermore,  even  in  equilibrium  where  fluid  motion  is  absent,  the  presence  of  a 
density  difference  between  the  solid  and  liquid  phases  can  have  thermodynamic  consequences. 

Typical  phase-field  models  for  a pure  material  without  convection  [7,  8,  9,  10]  introduce  a 
phase-field  variable  t),  which  is  a function  of  position  x and  time  t,  to  distinguish  the  solid 
from  the  liquid.  This  is  in  contrast  to  diffuse  interface  descriptions  of  a fluid  near  its  critical 
point  where  density  is  a more  natural  order  parameter  [11,  12].  The  phase-field  governing 
equations  have  been  derived  in  a thermodynamically-consistent  way  by  Penrose  & Fife  and 
others  [10,  13,  14]  following  the  formalism  of  irreversible  thermodynamics  [15].  AMW  followed 
these  same  principles  to  derive  their  phase-field  model  that  includes  convection.  They  retained 
the  use  of  the  phase-field  variable  <\>  and  allowed  the  density  of  each  phase  to  be  specified 
independently.  The  model  was  developed  to  describe  fluid  motion  that  is  incompressible  in 
the  bulk  phases.  Therefore,  in  order  to  account  for  expansion  or  contraction  flows  that  are 
driven  from  solidification  at  the  solid-liquid  interface  they  incorporated  the  ideas  of  quasi- 
incompressibility first  used  in  a diffuse-interface  description  of  a binary  fluid  by  Lowengrub  and 
Truskinovsky  [16]. 

AMW  treated  both  phases  as  fluids,  although  they  examined  the  case  where  one  phase  has  a 
sufficiently  large  viscosity  so  that  it  may  be  interpreted  as  an  approximation  to  a solid  phase.  In 
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many  solidification  applications,  a fluid  model  is  used  for  the  thermodynamic  description  of  the 
solid  phase,  in  that  the  elastic  properties  of  the  solid  are  ignored.  AMW  also  assumed  that  the 
phase  transition  was  first-order  and  had  an  anisotropic  surface  energy,  which  is  unconventional 
for  a fluid-fluid  system,  but  is  consistent  with  our  primary  intention  of  modeling  a solid-liquid 
system. 

There  have  been  a number  of  other  phase-field  type  descriptions  of  solidification  that  include 
convection  in  the  melt  (Caginalp  and  Jones  [17,  18],  Diepers  et  al.  [19],  Tonhardt  and  Amberg 
[20,  21]).  The  AMW  model  (see  also  reference  [22])  incorporates  in  a thermodynamically- 
consistent  way  the  density  effect  in  a phase-field  description  of  solidification  as  well  as  the 
appropriate  form  of  the  Korteweg  stress  term  in  the  momentum  equation. 

The  object  of  this  paper  is  to  examine  in  detail  the  governing  equations  for  the  phase-field 
convection  model  of  AMW  (see  also  reference  [22])  in  the  sharp-interface  limit  to  derive  the 
interfacial  conditions  of  the  associated  free-boundary  problem.  The  importance  of  this  analysis 
is  that  it  reveals  the  underlying  physical  mechanisms  built  into  the  phase-field  model  in  the 
context  of  a free-boundary  problem  and,  in  turn,  provides  a further  validation  of  the  model. 
To  date  sharp-interface  analyses  of  phase-field  models  have  concentrated  on  situations  in  which 
flow  is  not  present  and  the  densities  of  both  phases  are  the  same  [23]— [30] . The  version  of  the 
AMW  model  without  fluid  flow  has  recently  been  examined  in  detail  [30].  Therefore  the  focus 
in  the  present  work  is  on  the  effects  associated  with  density  differences  between  the  bulk  phases 
as  well  as  fluid  flow  and  how  they  manifest  themselves  in  the  sharp-interface  limit.  We  follow 
the  so-called  “classical”  sharp-interface  analysis  and  refer  the  reader  to  Fife  and  Penrose  [26], 
Karma  and  Rappel  [27,  28],  Almgren  [29]  and  McFadden,  Wheeler  and  Anderson  [30]  for  other 
possible  distinguished  limits  and  related  issues. 

The  paper  is  organized  as  follows.  The  governing  equations  are  described  briefly  in  Sec- 
tion 2 and  their  dimensionless  forms  used  in  the  present  analysis  are  given  in  Section  3.  The 
sharp-interface  analysis  is  first  examined  for  the  equilibrium  case  in  Section  4.  Section  5 dis- 
cusses scaling  issues  associated  with  the  nonequilibrium  version  of  the  equations  which  are  then 


3 


analyzed  in  two  different  limits.  Conclusions  are  given  in  Section  6. 


2.  Governing  Equations 


In  the  absence  of  convection,  phase-field  models  can  be  derived  by  extending  the  formalism  of 
classical  irreversible  thermodynamics  [15]  to  nonclassical  systems  whose  thermodynamic  poten- 
tials involve  square-gradient  energy  terms  [10,  13,  14].  AMW  built  on  these  ideas  to  develop  a 
phase-field  model  that  incorporates  convection.  The  governing  equations  of  this  model  which 
describe  convection,  heat  transfer  and  solidification  of  a pure  material  are 
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where  p is  the  density,  u is  the  velocity,  m is  the  stress  tensor,  </>  is  the  phase-field  variable  (we 
adopt  the  convention  that  cf)  — 0 denotes  the  liquid  phase  and  4>  = 1 denotes  the  solid  phase),  T 
is  the  temperature,  and  mg  : Vu  is  the  contraction  of  the  tensors  mg  and  Vu.  Additionally,  M 
is  a constant  that  represents  the  mobility  of  the  interface,  ep  is  the  Helmholtz  gradient  energy 
coefficient  given  by  c2F(T)  = e%  + Te|  where  e#  and  eg  are  the  gradient  energy  and  entropy 
coefficients  which  are  assumed  to  be  constants.  Further,  k((f>)  is  the  thermal  conductivity,  and 
g(T,  p,  <fi)  = e — Ts  + p/p  is  the  Gibbs  free  energy  per  unit  mass,  where  e is  the  internal  energy 
per  unit  mass,  s is  the  entropy  per  unit  mass  and  p is  the  pressure.  The  quantities  m and  mg 
are  given  by 
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where  r is  the  viscous  stress  tensor 


r = p{(ft) 


Wu+{Wu)T  - |(V-w)  / 


and  yu(</>)  is  the  viscosity  which  is  a function  of  (ft  given  by 


(3) 


p{<t>)  = Vsr{(ft)  + pL[  1 - r{(ft)]. 


(4) 


where  ps  and  pl  are  the  viscosities  in  the  bulk  solid  and  liquid  phases,  respectively.  The  function 
r(<ft)  is  monotonic  increasing  with  r(0)  = 0 and  r(l)  = 1;  suitable  choices  include  r(cft)  — (ft  or 
r((ft)  = (ft2(3  — 2(ft).  Here  T(V</>)  is  a homogeneous  function  of  degree  unity  that  permits  a general 
anisotropic  surface  energy  of  the  solid-liquid  interface,  and  £ is  the  generalised  ^-vector  whose 
components  are  defined  by  £j  = <9r(V (ft) / d(ftj,  see  Wheeler  and  McFadden  [31,  32].  We  will 
henceforth  confine  our  attention  to  case  of  isotropic  surface  energies  and  set  T(V</>)  = \W(ft\,  in 
which  case  £ = Weft/ |V0|. 

The  density  of  the  two  bulk  phases  may  be  different,  and  we  shall  assume  that  p depends  on 
(ft  alone, 

p{<t>)  = Psr{(ft)  + Pl[1  ~ r{(ft)},  (5) 


where  ps  and  pL  are  the  densities  in  the  bulk  solid  and  liquid  phases,  respectively.  This  as- 
sumption, in  which  the  density  does  not  depend  on  temperature  or  pressure,  known  as  quasi- 
incompressibility, places  a constraint  on  the  form  of  the  underlying  thermodynamic  potentials 
[16]  and  requires  that  the  underlying  Gibbs  free  energy  density  be  given  by 

9(T,p,0)  = 9„(r»+^T^'  <6) 

w) 

where  pr  is  a reference  pressure  associated  with  an  isothermal  stationary  planar  interface  at 
melting  temperature  Tm-  We  take  go{T,  (ft)  to  have  the  form 
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in  which  case  the  corresponding  expressions  for  the  internal  energy  and  entropy  densities  are 
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where  1 /a#  = 1/a  — 1/as ■ Here  1/a  is  the  height  of  the  double  well  of  the  Gibbs  free  energy 
density  at  T — Tm , and  1/a#  and  1/as  are  the  heights  of  the  double  wells  in  the  internal  energy 
and  entropy  densities,  respectively.  The  quantity  eo  is  a constant  reference  energy  and  both  the 
heat  capacity  per  unit  mass  c and  the  latent  heat  per  unit  mass  L are  assumed  to  be  constant. 

The  double  well  function  Hm(4 >)  appears  in  the  “per  unit  mass”  quantity  go(T,  (f)).  The 
analog  “per  unit  volume”  quantity  associated  with  p(4>)go{T,4>)  is  p{(j>)Hm{<t))  = Hv{<j>).  A 
common  form  of  a double-well  potential  is  the  function  4>2{1  — 4>)2,  and  in  standard  phase-field 
models  that  do  not  include  convection  it  is  usually  associated  with  a “per  unit  volume”  quantity 
(e.g.  Wang  et  al.  [10]).  In  applications  in  which  the  density  of  the  two  phases  are  constant  and 
equal,  the  per  unit  mass  and  the  per  unit  volume  specification  of  the  double  well  are  equivalent. 
However,  in  the  present  situation  the  bulk  densities  are  not  necessarily  equal  and  therefore  we 
shall  examine  the  merits  of  the  per  unit  mass  form  Hm{(j))  = 4>2{  1 — 4>)2  and  the  per  unit  volume 
form  Hv(<t>)  = Pl4>2{  1 - 4>)2  [in  which  case  Hm{<j>)  = Pl4>2{ 1 - <i>)2  / p(<i>)]. 

The  one  dimensional  solution  of  the  governing  equations  for  an  isothermal  system  at  the 
melting  temperature  Tm  with  equal  solid  and  liquid  densities  (ps  = Pl)  is  given  by 


4> 


1 
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1 — tanh 
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3to  , 4 / 2 

p = Pr ^rsech  — 

y 161  \2l 


(10) 


The  interface  thickness,  /,  and  surface  energy  70,  are  related  to  the  phase-field  parameters  by 

(11) 


, ,rri  \ /2a  Cf(Tm)  fpL 
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It  is  also  convenient  to  define  an  associated  capillary  length  by  lc  = TmIq/ {Pl[L2 / c]). 


3.  Dimensionless  Governing  Equations 

We  non-dimensionalise  the  governing  equations  by  introducing  a reference  length  scale  Iq,  which 
is  a typical  length  scale  associated  with  the  interface  shape,  such  as  a dendrite  tip  radius.  The 
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reference  scales  for  time,  velocity,  pressure  and  stress  are  based  on  thermal  diffusion  in  the 
liquid.  The  temperature  scale  is  L/c  and  the  density,  viscosity  and  thermal  conductivity  are 
scaled  with  their  values  in  the  bulk  liquid.  The  dimensionless  quantities,  denoted  by  primes,  are 
related  to  their  dimensional  counterparts  by 
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where  k,l  is  the  thermal  diffusivity  in  the  bulk  liquid  phase.  The  dimensionless  governing 
equations  are 
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where,  for  simplicity,  we  have  omitted  the  primes  on  the  dimensionless  variables.  The  dimen- 
sionless internal  energy  density,  entropy  and  Gibbs  free  energy  are  given  by 
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where  e'0  = eo/L.  Note  that  g is  scaled  with  l/(2a),  s with  c/(2aL)  and  e — e0  with  L.  The 
dimensionless  tensors  are 


m — ap  + (l  + aO)^  + Prr, 
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The  source  term  in  the  energy  equation  is 
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The  dimensionless  parameters  are  given  by 
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where  27,  is  the  kinematic  viscosity  of  the  liquid  phase  and 


(24) 

(25) 

(26) 


W)  = 1 + (jr  - rW>  f*W  = 1 + ^(<^>),  p(4>)  = 1 + (^  - r(0).  (27) 

We  note  that  the  parameter  7 is  related  to  the  capillary  (or  crispation)  number  Ca  by  7 = 
6 Pr/Ca,  where  Ca  = hlk>l/{1  o7o)-  Additionally,  if  one  follows  the  thermodynamic  relationships 
between  the  gradient  coefficients  and  the  double- well  heights,  namely,  c2f(Tm ) = e2E  + Tm£2s  and 
1/oe  = 1/a  — 1/as,  we  find,  as  in  McFadden,  Wheeler  and  Anderson  [30]  that  a , (3,  S and  v are 
not  independent  parameters  but  rather  are  related  through 


/3  — a = A (u  - 5), 
a — S — \v. 


(28a) 

(28b) 


In  the  absence  of  flow  these  equations  reduce  to  the  generalised  phase-field  equations  recently 
studied  in  Refs.  [30]  and  [33].  We  observe  that  e is  the  ratio  of  the  interface  thickness  l to  the 
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reference  length  scale  l0  and  in  the  analysis  that  follows  we  shall  consider  the  sharp-interface 
limit  e — > 0.  The  leading-order  free-boundary  problem  that  emerges  from  a sharp-interface  limit 
of  these  equations  depends  on  the  distinguished  limit  that  is  taken  [23,  27,  30].  The  particular 
sharp-interface  limit  in  which  A = 0(1)  as  e — » 0,  where  the  interface  thickness  is  comparable 
to  the  capillary  length,  is  the  so-called  ‘thin  interface’  limit  studied  by  Karma  and  Rappell  [27]. 
The  case  A = O(e)  as  e — * 0,  where  the  interface  thickness  is  much  smaller  than  the  capillary 
length,  is  the  so-called  ‘classical’  limit. 

In  this  paper  we  consider  a classical  sharp- interface  limit  for  the  simplified  model  where 
a — (3  = 0 and  v = 5 = S/ A.  Here  the  square-gradient  and  double  well  terms  are  present  in  the 
internal  energy  functional  but  not  in  the  entropy  functional.  Notice  that  since  A = 0(e)  in  the 
classical  limit,  we  have  v — S — 5/A  = 0(e~ !);  these  terms  all  appear  in  the  energy  equation. 
We  refer  the  reader  to  McFadden,  Wheeler  and  Anderson  [30]  where  the  more  general  case  of 
nonzero  a,  (3,  5 and  v is  examined  without  fluid  flow.  The  dimensionless  governing  equations  in 
this  case  are 

^ + V • (pu)  = 0,  (29) 
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(32) 


where  m and  H are  given  by  equations  (21a)  and  (23)  with  a = 0. 


4.  Equilibrium:  Sharp-Interface  Analysis 

In  this  section  we  conduct  a classical  sharp-interface  analysis  of  the  equilibrium  version  of 
equations  (29)-(32).  Here  the  fluid  velocity  u and  all  time  derivatives  are  zero.  Additionally, 
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the  energy  equation  admits  an  isothermal  solution  where  the  temperature  is  denoted  by  6.  In 
the  classical  limit,  e — > 0 with  A = eA  and  A = 0(1),  thin  interfacial  regions  of  thickness  O(e) 
develop  in  cf).  Away  from  these  regions  the  system  is  in  the  solid  or  liquid  bulk  phase  and  we 
express  the  solution  in  these  so-called  outer  regions  as  regular  perturbation  series  in  powers  of 
e.  We  write 


p(C\C2,C3) 

— Po(C\C2,  (3)  + ePi((\  Ci  C3)  + • • • 5 

(33a) 

0(c\c2,c3) 

— 0o(C\  C2>  C3)  + ^(C1,  C2,  C3)  + • • • 5 

(33b) 

where  C2?  C3  are  coordinates  locally  body-fitted  to  the  mean  position  of  the  interface,  given  by 
the  surface  <j)  = 1/2.  Here  £2  are  surface  coordinates  measuring  distance  along  the  principal 
directions,  and  £3  is  a coordinate  measuring  distance  in  the  direction  normal  to  the  surface 
4>  — 1/2,  which  may  also  be  expressed  by  = 0.  We  adopt  the  convention  that  C3  is  positive  in 
the  liquid  and  negative  in  the  solid.  Further  details  of  the  differential  geometry  of  body-fitted 
coordinates  are  given  in  Appendix  A. 

The  phase-field  equation  in  the  outer  regions  gives  that  0O  = 0 in  the  liquid  and  </>0  = 1 in 
the  solid.  The  momentum  equation  V • m = 0 in  the  outer  region  requires  that  Po(Cl,  C, ’2,  C3)  is 
constant  in  each  phase.  We  put 

Po  = Ps  ~ P*->  Po  = Pl  ~ P*i  (34) 

where  p*s  = PsIq/ {Pl^l),  Pl  — Pl^q/ {Plk2l),  and  ps  and  pi  are  the  dimensional  pressures  in  the 
bulk  solid  and  liquid  phases,  respectively. 

We  now  consider  the  solution  in  the  interfacial  layer.  We  observe  from  the  one-dimensional 
solution  (10)  that  the  dimensional  pressure  in  the  interface  is  characterized  by  q0 //,  and  hence 
the  dimensionless  pressure  in  the  interface  scales  as  (to/OA/^lAo)  — 7/(6e)  where  7 = 0(1). 
Consequently,  we  write  the  inner  expansions  as 

PiCX2,0  = )[Po«I.C2,C)  + ^i(C1,C2.0  + .-- 
<KC’,C2.  C)  = $o(C\<2,0  + e'MC'.CsO  + ■ • • ; 


10 


(35a) 

(35b) 


here  and  below  we  use  capital  letters  to  denote  the  inner  variables  and  ( = £3/e  is  the  normal 
coordinate  scaled  to  the  thickness  of  the  interfacial  layer.  Equilibrium  conditions  follow  from 
the  momentum  and  phase-field  equations  as  shown  below. 


4.1  Momentum  Equation 


We  examine  the  first  two  orders  of  the  momentum  equation  0(e  2)  and  0(1)  to  find  that 

2 


Po  = c0(C\C2)  ~ t:7 


1 / d$c 


1 X2\  „*  , \ /c  (d*oY  u 


d( 


2 ' v d(  J ’ 

Pi  — ci(C\ O - p*  + 7(^1  + ^2)  [ ( 

*/— 00  y 

where  c0  and  Ci  are  independent  of  £. 

Matching  the  inner  and  outer  pressures  requires 

0=  lim  jp0(r0)  - -(P0  + cPi)  + . . .)  , 

C->±oo  L e J 

and  gives  co(£\C2)  = 0 and  cdC1,^2)  = p*s.  Additionally, 


d(  - 7- 


5C  <9£ 


Po 


(36a) 

(36b) 


(37) 


(38) 


which  relates  the  jump  in  the  outer  pressure  across  the  interface  to  the  interface  curvature.  We 
identify  the  integral  factor  below  after  treating  the  phase-field  equation. 


4.2  Phase-Field  Equation 

The  0(1)  inner  problem  for  the  phase-field  equation  can  be  written  as 


0 = y,ll°  | 1 d 


ln(p(*o)) 


(39) 


dC2  2 d<\> 

where  we  have  made  use  of  the  expression  (36a)  for  P0.  This  equation  can  be  simplified  if  we 
first  multiply  through  by  the  integrating  factor  (l/p($ 0))1//2  to  obtain 

( 1 Y p(^)H'm(%) 


( 1 

U$o)/  ac 


,p($o) 


(40) 


11 


and  then  multiply  this  result  by  (l/p($o))^2^o/^C-  Integration  then  gives  the  first  integral 

2 


0 = 


d$0 

dC 


p($o)*U*o). 


(41) 


The  0(e)  inner  phase- field  problem  can  be  expressed  as 


d 2$1  d 
dC  + d<\> 


In 


1 


A®  o) 


<9$0  d$i 


In 


A*  0)) 


Po 

7 


= Ru 


(42) 


where 


Ri  = (1C,  + £2)^  + A9/>($„)r'($o)  + - |r 

<7(,  7 <70 


In 


p($o) 


+ (/Cl+C2)^ 


In 


1 


p($o) 


r.r 


<9C 


dC, 


(43) 


and  we  have  made  use  of  equation  (36b)  for  Pi.  Equation  (42)  may  be  manipulated  to  obtain 
the  solvability  condition 

r+00  1 <9$0 


0 


= /_' 


■Pi^C, 


-00  p($o)  d( 

which  can  be  simplified  to  obtain  the  Clausius-Clapeyron  relation 


(44) 


0 = -A6  + 


(Ps  ~P*)  (1  PL 


t1  ~ i] + {Ki + ^ L (l2) dc 


(45) 


7 V Ps  J J~ 00  V 

This  equation  expresses  temperature  in  terms  of  the  bulk  pressure  and  interface  curvature. 

In  dimensional  form,  the  pressure  jump  condition  (38)  and  the  curvature-modified  Clausius- 
Clapeyron  relation  (45)  are 


Ps  = 7(^1  +£2), 

0 = L ( 1 — + ( ps  - Pr)  f \ + — (/Cx  + /C2), 

V TmJ  \pi  ps ) Pl 

where  7 is  the  surface  tension  (or  excess  Kramers  potential)  defined  in  Appendix  B as 

7 = 670  f ^ d^  = 6l0J0 


(46a) 

(46b) 


(47) 
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and  is  determined  by  equation  (41).  Equations  (46)  determine  the  liquid  pressure  pL  and 
the  temperature  of  the  system  T in  terms  of  the  solid  pressure  p5,  the  interface  curvature,  and 
the  material  parameters. 

The  value  of  7 is  not  in  general  equal  to  70  defined  by  equation  (11)  but  rather  depends 
on  the  density  ratio  ( ps/pL ) and  the  form  of  p(0)i/m(0).  When  p(0)ifm(0 ) = 4>2{  1 - </>)2  [so 
that  the  double-well  potential  is  associated  with  a “per  unit  volume”  quantity,  see  equation  (7)] 
equation  (41)  admits  the  standard  hyperbolic  tangent  profile  and  the  integral  in  equation  (47) 
evaluates  to  1/6  so  that  7 = 70 . However,  when  Hm{4>)  = 02(1  — (f>)2  [so  that  the  double- well 
potential  is  associated  with  a “per  unit  mass”  quantity,  see  equation  (7)]  equation  (41)  no  longer 
has  the  hyperbolic  tangent  profile  solution.  Additionally,  with  r(<f>)  = </>2( 3 — 2 (j>)  we  find 


(48) 


and  with  r(0)  — <\>  we  find 


7 ; 7o- 


3 + 9 (^)  ~2  + 11^  + 9 (^Y  +3  (^Y 

\PL  I PL  \PL  J \PL  ) 


35 


1 + 


(ft) 


(49) 


Each  of  the  expressions  for  7 in  (48)  and  (49)  evaluates  to  70  when  the  bulk  densities  are  equal. 
Furthermore  they  are  increasing  functions  of  Ps/pl  and  vary  approximately  from  O.770  to  2.37o 
as  ps/pl  varies  from  0.1  to  10. 

This  analysis  shows  that  the  equilibrium  version  of  the  present  phase-field  model  allows  an 
isothermal  system  in  which  the  pressure  difference  between  the  two  bulk  phases  is  equal  to  the 
surface  tension  7 times  interface  curvature  (the  Young-Laplace  equation)  while  the  temperature 
of  the  system  is  determined  by  7 times  interface  curvature  (the  Gibbs-Thomson  effect)  and  also 
by  the  bulk  pressure,  when  there  is  a density  mismatch,  through  a modified  form  of  the  Clausius- 
Clapeyron  relation.  The  Clausius-Clapeyron  pressure  effect  has  been  investigated  by  LaCombe 
et  al.  [34]  as  an  additional  control  on  the  interface  melting  temperature  during  dendritic  growth. 
We  also  note  that  Maruyama  et  al.  [35]  have  examined  transitions  in  the  kinetic  growth  shapes 
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of  ice  7/j,  from  a circular  disk  to  a hexagonal  plate,  in  response  to  pressure-induced  alterations 
of  the  melting  point  near  the  roughening  transition. 


5.  Nonequilibrium:  Sharp-Interface  Analysis 

In  this  section  we  examine  the  role  of  nonequilibrium  effects  on  the  hydrodynamic  boundary 
conditions  that  arise  in  the  sharp-interface  limit  as  well  as  the  heat  flux  and  Clausius-Clapeyron 
boundary  conditions. 

5.1  Scaling  Arguments 

In  order  to  guide  this  analysis,  we  examine  the  Clausius-Clapeyron  relation  for  a planar  interface 
derived  in  AMW  [1],  In  terms  of  the  present  nondimensionalizations,  their  result  is 


-fA e-e-  ^-i  (pj-p*) 
7 PS 


The  left-hand  side  of  this  equation  contains  the  temperature  and  pressure  terms.  The  right-hand 


side  involves  an  interfacial  attachment  kinetics  term,  a volume-change  term  (which  can  also  be 


interpreted  as  a velocity-dependent  modification  of  the  far  field  pressure  owing  to  the  vapor 
recoil  effect  [36])  and  a viscous  dissipation  (or  viscous  drag)  term,  respectively.  The  attachment 


kinetics  term  is  present  in  phase- field  models  without  fluid  flow  [10]  although  the  integral  term 


here  is  modified  by  the  density  function.  The  viscous  drag  term  is  associated  with  the  velocity 
change  across  the  interfacial  region  in  the  normal  direction  and  vanishes  for  inviscid  flows  or 
density-matched  phases.  In  the  analysis  of  solute  trapping  for  a binary  alloy  by  Ahmad  et 
al.  [37],  an  analogous  effect,  called  solute  drag,  was  identified  [see  the  integral  term  in  their 
equation  (34)].  In  equation  (50)  the  viscous  drag  term  vanishes  with  the  square  of  the  density 
difference  across  the  interface  times  the  interface  velocity  while  the  solute  drag  effect  in  Ahmad 
et  al.  vanishes  with  the  square  of  the  concentration  difference  across  the  interface  times  the 
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interface  velocity. 

In  an  asymptotic  sense,  it  is  possible  to  obtain  a number  of  dominant  balances  among  the 
physical  effects  represented  by  equation  (50).  In  the  case  of  equal  densities  and  no  flow  possible 
balances  include  the  situations  where  temperature  is  controlled  at  leading  order  by  attachment 
kinetics  (e.g.  if  the  interface  is  nearly  planar),  by  curvature  effects  (the  Gibbs  Thomson  effect), 
by  both  effects,  or  by  neither  (e.g.  see  references  [10,  23,  27,  28]).  In  the  present  context,  we 
must  determine  how  density  differences  and  fluid  flow  affect  the  interface  temperature.  If  in 
equation  (50)  all  parameters  are  0(1)  as  e — » 0,  then  the  viscous  drag  term  dominates  and  no 
balance  can  be  achieved.  Hence  we  now  reconsider  the  different  terms  in  the  equation  in  more 
detail. 

Anderson,  McFadden  and  Wheeler  [1]  found  that  for  a typical  material  such  as  lead  the 
viscous  drag  term  was  negligible  compared  to  the  attachment  kinetics  term  due  at  least  in  part 
to  a small  density  difference.  This  was  true  even  for  solid  to  liquid  viscosity  ratios  of  as  large 
as  106.  Consequently,  in  Section  5.2  we  consider  an  asymptotic  regime  in  which  the  density 
difference  ( Ps/Pl  — 1)  — * 0 as  e — » 0. 

In  Section  5.3  we  consider  a second  asymptotic  regime  in  which  the  density  difference  is  0(1) 
and  the  coefficient  Pr/y  = Ca/6  = O(c)  while  the  viscosity  ratio  ps/ is  treated  as  0(1)  with 
respect  to  e.  This  scaling  may  be  more  appropriate  for  two  fluids  (say  liquid  and  vapor)  where 
the  density  difference  is  not  small  and  viscous  effects  are  of  secondary  importance. 

Another  possible  scaling  would  be  to  specify  the  viscosity  function  p(cf>)  to  be  0(e)  in  the 
interfacial  region  where  dp/dC,  is  nonzero.  However,  it  is  not  clear  that  this  is  physically  plausible 
and  we  do  not  pursue  this  further. 

As  a simplification,  in  the  remaining  analyses  we  shall  use  exclusively  the  double-well  po- 
tential given  by  Hv(4>)  = p((f))Hm((j))  = cf)2(l  — <j))2.  We  note,  however,  that  in  the  case  of  nearly 
matched  densities  in  Section  5.2,  Hv(<p)  and  Pm(</>)  are  equivalent  at  leading  order. 
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5.2  The  Case:  ps/pL  ~ 1 = O(e),  Ca  = 0(1),  A - O(e) 

In  this  section  we  describe  a classical  sharp-interface  analysis  where  A = eA  with  A = 0(1)  as 
e — > 0 and  we  assume  that  the  solid  and  liquid  densities  are  nearly  matched.  In  particular,  we 
take 


p(<j>)  = l+exr(<f>),  (51) 

where  x is  0(1)  as  e — » 0. 

We  use  the  expansions 

u(C\C2,<3,«)  = «o(C\  C2,  C3,  <)  + «5i(c‘ , c2,  c3.  <)  + ••• , (52a) 

P(C\C2,C3,«)  = Po(C1,C2.C3.*)  + epi(C1X2,<3.«)  + ---.  (52b) 

0(C'.C2,C  3,<)  = MC\C2X3,t)+eMCl,(2X3,  <)  + ...,  (52c) 

0(C\C2,C  3,t)  = e0(C1,C2.C3,<)  + e«i(C1.<2.C3.t)  + .--.  (52d) 


for  the  outer  variables.  The  leading  order  velocity  and  pressure  in  the  bulk  phases  satisfy  the 
incompressible  Navier-Stokes  equations 


Duo 

Dt 

0 


— Vp0  + TrV  • r0, 


V • u0. 


(53a) 

(53b) 


The  phase-field  equation  gives  <po  = 0 in  the  liquid  and  4>q  = 1 in  the  solid  and  the  energy 
equation  becomes 


D60 

Dt 


QPr 

= fcV2f>0  + -p-ro  : V«o, 

Ay 


(54) 


where  k = 1 in  the  liquid  and  k = ks/ki  in  the  solid. 

The  outer  variables  must  satisfy  interfacial  conditions  which  we  will  determine  by  examining 
the  inner  problem  and  invoking  matching  and  solvability  conditions.  In  the  inner  region  we  put 


£?(C,C2,C,*)  = Vo(C1.C2-C,<)  + e£?i(C1,C2,C,0  + ---. 


(55a) 
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P(C\C2,C,i)  = J [foK\C2,C,<)  + eP1(C1,C2,C, ()  + •••]  , (55b) 

m sc2,  c,*)  = 4>o(c1,c2.c*)  + f<Mc1,c2,c, (55c) 

®(C\C  2iCit)  — 0o(C\  C2i  C»  £)  + €@i(C1>  C2>  Ci  t)  + • • • , (55(1) 

M(C\C2,C,i)  = -e  + £M1«1,C2,C, ()  + ...].  (55e) 

We  examine  each  of  the  governing  equations  in  the  following  sections. 


5.2.1  Continuity  Equation 

The  leading  order  problem  for  the  continuity  equation  (29)  shows  that  Uq  is  independent  of  (. 
The  matching  condition  for  velocity  is 


lim  Uq  — m0|l,  lim  U0  = u0\s,  (56) 

£— >■  00  £ — ► — 00 

where  ‘|  p and  ‘Is’  denote  a quantity  in  the  bulk  phase  evaluated  on  the  liquid  and  solid  side  of 
the  interface,  respectively.  It  follows  that 


u0  • n 


= 0. 


(57) 


This  is  the  normal  mass  flux  condition  when  the  densities  match  to  leading  order. 


5.2.2  Momentum  Equation 

At  leading  order  [C9(e-1)]  the  momentum  equation  (30)  gives  that 

dM31  _ OMq2  _ <9Mg3 
d(  d(  ~ d(  ~ ’ 

Using  (21a),  (22)  and  the  results  in  Appendix  A we  find  that 


(58) 


M031  = Prr“\  M32  = Pr  rf,  = - 


31  A//32  _ ~-32  71  /T 33  _ 


1 _ / d<S>0 

F°  + 27  l 9C" 


+ Pr  r33. 


(59) 


Integrating  (58)  once  gives  that  these  leading  order  stresses  are  independent  of  £.  However, 
applying  the  matching  condition  on  the  stress  gives 


lim  M o = 0 

C— >±oo 


(60) 
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and  so  Mg 1 = Mg2  = Mg3  = 0.  Employing  the  expressions  for  Tg1,  t32  and  Tq3  given  in  Appendix 
A we  find 


»(«f  -o,  .0.  - 


p , 1 ~ (d®° 

Po  + 27  ( sF 


+ = 0, 


(61) 


so  that  the  forms  of  the  velocity  components  and  the  pressure  in  the  interfacial  layer  are 


U'a  = MC‘,C2,*),  y?  = MCI.C2a).  C?  = MCtC2.*). 


p„ = -57 


9C 


(62) 

(63) 


where  62  and  63  are  independent  of  £.  The  matching  condition  for  velocity  (56)  then  gives 
the  no-slip  condition 


[So  - "(So  - ft)]  s = 0, 


(64) 


which  shows  that  the  tangential  components  of  the  leading-order  velocities  in  the  bulk  phases 
are  continuous  across  the  interface. 

Consistent  with  the  leading-order  inner  stresses  vanishing  in  the  far-field  limit,  the  leading- 
order  inner  pressure  P0  also  vanishes  in  this  limit.  We  emphasize  that  although  the  leading 
order  inner  pressure  term  is  0(1 /e)  [see  equation  (55b)],  the  leading  order  contribution  to  the 
stress  is  0(1)  as  M0  = 0.  Stress  jumps  across  the  interface  are  determined  next. 

The  momentum  equation  at  0(1)  in  the  inner  region  gives 


(U30-Vn)^- 

M-v.)  f 


- (Ki  + KyM„31  - /Ci(M013  + M31),  (65a) 

dMf2  dMf2  dM22 
~dC~  + ~dC~+  dC2 
dM33  dMf3  dM23 


d(  dc 


- (/C,  + JC2)M32  - /C]_(Mq3  + M32),  (65b) 
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r32' 


(/C1  + /C2)M33  + A:1Mo11  + /C2M22.  (65c) 
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r22 


dC 


However,  using  Mg1  = M32  = Mf3  = 0,  the  symmetric  nature  of  the  stress  tensor,  equa- 
tions (62),  and  the  expressions  for  the  components  of  the  stress  tensor  given  in  Appendix  A, 
these  equations  simplify  and  we  find  that 


dM31 


d_ 

dC 


Po + 27 


(66a) 
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(66b) 


dMf 2 

~dT 

dMf 3 

~dC~ 


_d_ 

dc2 


p , ( d$o 

~p°  + 27  ( ' ~ac 


= -(Kj+Kj) 


-Fo  + 27by 


(66c) 


Integrating  these  equations  from  ( = — oo  to  £ = -foo  gives  that 


M,31'  M? 

-oo 


+ 0O  _ <9X  . ,32  +00  _ dX  . .00  +00 


= -7a?' 


33 


= -7(a:,  + k2)x, 


(67) 


where,  on  using  equation  (63)  the  integral  X may  be  expressed  as 


(68) 


The  integral  X evaluates  to  1/6  for  the  double-well  Hv{(j>)  = p((f))Hm((f))  = </>2(l  — 8>)2.  Matching 
allows  equations  (67)  to  be  expressed  in  dimensional  form  as 


[{-Pol  + T0)  • n]x 


[to  • n\ || 


L 

5 

L 

S 


— “7o(^l  + £2), 

= 0, 


(69a) 

(69b) 


where  the  subscripts  ± and  ||  indicate  components  normal  and  parallel  to  the  interface.  Here 
To  represents  the  leading-order  outer  viscous  stress  tensor.  These  are  the  conventional  stress 
balances  that  hold  in  the  absence  of  Marangoni  effects  whereby  the  tangential  stress  is  continuous 
across  the  interface  and  the  normal  stress  difference  across  the  interface  is  balanced  by  surface 
tension  times  interface  curvature. 

In  the  interest  of  modeling  a solid-liquid  system  we  examine  the  limit  ms/Ml  — > 00.  In  the 
bulk  solid  this  limit  gives  at  leading  order  that  V • (Vus  + Vu^)  = 0 while  the  above  stress 
boundary  conditions  give  [(Vus  + Vi/s)  • n]||  = 0 and  [(Vus  + Vu^)  • n]x  = 0.  Hence  the  solid 
phase  undergoes  rigid  body  motion  at  leading  order.  The  leading  order  velocity  field  in  the 
liquid  phase  satisfies  the  Navier-Stokes  equations  in  the  bulk  and  is  subject  to  the  normal  flux 
and  no  slip  boundary  conditions  on  the  solid-liquid  interface.  As  expected  in  this  limit,  the 
stress  in  the  liquid  is  not  imposed  through  a boundary  condition  but  rather  can  be  derived  from 
the  velocity  field. 


19 


5.2.3  Energy  Equation 


We  next  consider  the  energy  equation  (32)  and  recall  that  v = 6 = (5/e A)  = We  note 

that  e = 6 - r(</>)  + (S/2Ae)Hm((j))  — ( Sp*/Ayp ),  so  that  in  the  inner  region  the  advective  term 
in  the  energy  equation  contributes  at  leading-order  0(e~2).  The  effects  of  dissipation  as  well  as 
thermal  diffusion  also  contribute  at  this  order  and  we  find  that 
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We  note  that  the  leading  order  phase- field  variable  <f>0  satisfies  equation  (74)  so  that  terms  with 
coefficients  5/2A  cancel.  Equation  (70)  can  then  be  integrated  twice  to  give 
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where  d0  and  d\  are  independent  of  £.  Matching  the  inner  temperature  to  the  outer  solutions 
gives  d\  — 0 so  that  0o  = d0(C\  C2>  t)  and  0o\l  = #o|s  = 0o-  Therefore  the  inner  temperature  0o 
is  independent  of  C,  and  the  outer  temperature  is  continuous  across  the  interface.  In  dimensional 
terms  we  write  Tl  = Ts  = Tj  where  Tj  is  the  interface  temperature  as  given  in  equation  (78). 

We  can  identify  a heat  flux  boundary  condition  when  we  examine  the  0(e~1)  problem.  The 
analysis  involves  the  inner  velocity  component  correction  C/f , which  is  determined  from  the 
correction  to  the  continuity  equation.  The  details  of  this  calculation  are  contained  in  Appendix 
C.  We  find  that  in  dimensional  form  the  heat  flux  boundary  condition  obtained  is 
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This  is  the  classical  Stefan  condition  modified  to  account  for  motion  in  the  solid  and  liquid 
bulk  phases  as  well  as  the  curvature  of  the  interface.  The  effect  of  flow  only  enters  in  the  first 
factor  where  it  arises  as  the  normal  velocity  of  the  material  relative  to  the  interface.  The  cur- 
vature effect  follows  from  the  term  DtW2(f)  in  equation  (32).  It  represents  the  internal 

energy  gradient  and  double-well  terms  in  our  model.  Similar  curvature  effects  have  been  recog- 
nised by  Umantsev  and  coworkers  [38,  39]  who  identified  the  coefficient  of  the  curvature  to  be 
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7o  — T d'ye / dT.  Equation  (72)  is  in  agreement  with  this  work  when  we  observe  from  (11)  that  70 
is  independent  of  temperature.  Fife  and  Penrose  [26]  and  Fried  and  Gurtin  [40]  have  also  iden- 
tified this  curvature  effect  in  the  context  of  phase-field  models.  In  sharp-interface  formulations 
Wollkind  and  Maurer  [41],  Umantsev  and  Davis  [38],  Zhang  and  Garimella  [42]  and  Schlitz  and 
Garimella  [43]  studied  its  effect  on  the  stability  of  the  interface  and  Lemieux  and  Kotliar  [44] 
its  influence  on  velocity  selection  during  dendritic  growth. 


5.2.4  Phase-Field  Equation 

At  leading-order,  the  phase- field  equation  (31)  for  nearly  matched  densities  becomes 
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This  equation  has  the  first  integral 
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The  first  order  phase-field  equation  is  given  by 
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where  we  have  used  (63)  and  (74).  Upon  multiplying  equation  (75)  by  <9<f>o/<9C  the  solvability 
condition  can  be  readily  obtained 
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With  Hm($ 0)  = $§(1  - <4>0)2  and  r(0)  = 0 and  r(l)  = 1 we  find  that  the  interfacial  condition  is 
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or  in  dimensional  terms 
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where  T/  is  the  interface  temperature  and  p,m0b  is  the  interface  mobility  given  by 
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Equation  (78)  gives  the  interface  temperature  Tj  in  terms  of  curvature  and  attachment 
kinetics  and  is  equivalent  (when  uo|s  — 0)  to  the  result  for  no  flow  [10].  It  differs  from  the  more 
general  Clausius-Clapeyron  relation  [see  equation  (102)]  in  that  the  effects  of  pressure,  density 
change  and  viscous  drag  do  not  enter  at  leading-order.  We  expect  that  this  equation  represents 
typical  solid-liquid  systems  in  which  the  density  difference  between  the  two  phases  is  small. 


5.3  The  Case:  Ca  --  O(e),  ps/ Pl  — 1 = 0(1),  A = 0(e) 


As  in  the  previous  case,  we  conduct  a classical  sharp-interface  analysis,  but  now  we  assume  that 
the  density  difference  is  0(1)  and  the  viscous  effects  are  O(c).  Specifically,  we  put  Pr/ 7 = 
Ca/6  = Ce  where  C = 0(1),  and  again  take  A = Ae  where  A = 0(1)  as  e — » 0. 

We  follow  the  same  procedure  as  in  the  previous  case  and  find  that  the  leading-order  outer 
problem  remains  the  same  with  the  exception  that  the  Navier-Stokes  equation  (53a)  is  replaced 
by  the  Euler  equation 
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where  in  this  equation  p — 1 in  the  liquid  and  p = ps/ Pl  in  the  solid  and  the  energy  equation 
is  without  the  dissipation  terms 
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where  k = 1 in  the  liquid  and  k — ki/ks  in  the  solid. 

Interfacial  conditions  are  derived  by  examining  the  inner  region.  Here  we  again  use  the 
expansions  (55)  and  consider  each  equation  in  turn. 
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5.3.1  Continuity  Equation 


The  leading  order  problem  for  the  continuity  equation  (29)  may  be  solved  to  obtain  the  form 
for  Uq  through  the  interface  as 
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The  matching  condition  (56)  applied  to  the  normal  velocity  component  gives 


p{u0  • h - Vn)  s = 0, 


(83) 


in  dimensional  form,  where  J0  = Pl{uq\l  • n — Vn)  — ps{uo\s  • n - Vn)  is  identified  as  the  mass 
flux  across  the  interface. 


5.3.2  Momentum  Equation 


At  leading  order  [0(e_1)j  the  momentum  equation  (30)  gives  equations  (58)  as  in  the  analysis 
of  Section  5.2.  However,  owing  to  the  different  scaling,  equations  (21a),  (22)  and  the  results  in 
Appendix  A now  lead  to 
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Equations  (58)  and  the  matching  condition  on  the  stress  (60)  require  Mq3  = 0 so  that 
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The  next  order  momentum  equation  in  the  inner  region  gives 
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However,  using  Mg 1 = Mg2  = Mg3  = 0,  the  symmetric  nature  of  the  stress  tensor  and  the 
expressions  in  Appendix  A with  Pr/7  = Ce  these  equations  simplify  to  give 
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Now  using  these  expressions  and  (85),  equations  (87)  become 
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The  leading-order  phase-field  d>0  satisfies  the  equilibrium  equation  (41)  [see  also  Section  5.3.4 
below]  and  depends  only  on  £.  Consequently,  we  integrate  equations  (89a)  and  (89b)  to  obtain 
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where  6X,  62,  63  and  64  are  independent  of  £.  We  now  apply  the  matching  condition  for  veloc- 
ity (56)  and  find  63  = 64  = 0 so  that  Uq  = 6i(C1,C2^)  and  Uq  ~ MCSC2^)-  Hence,  the  leading 
order  outer  tangential  velocity  components  are  continuous  across  the  interface 
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These  results  give  that  Mf1  = Mf2  = 0.  However,  integrating  the  momentum  equation  (89c) 

gives  that  Px  is 
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where  C3  is  independent  of  ("•  We  next  apply  the  matching  principle  for  the  pressure  (37)  to  find 

that  C3  = P0I5  + P*  + Jo  W0I5  • n and  the  jump  in  the  outer  pressure  at  the  interface  is  given  (in 
dimensional  form)  by 
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where  J0  is  given  in  the  text  after  equation  (83).  Therefore,  the  pressure  undergoes  a jump 
across  the  interface  due  to  the  vapor  recoil  effect  and  interface  curvature.  Viscous  stresses  do 
not  appear  at  this  order  in  the  analysis  since  Ca  = O(e). 

5.3.3  Energy  Equation 

We  examine  the  energy  equation  (32)  and  find,  as  in  the  previous  case,  that  TL  = Ts  = 77; 
the  outer  temperature  is  continuous  across  the  interface  with  the  interface  temperature  Tj  given 
by  equation  (102).  The  heat  flux  boundary  condition  follows  in  an  analogous  manner  as  that 
of  Section  5.2.3.  The  details  of  this  calculation  are  given  in  Appendix  D.  We  find  that  in 
dimensional  form  the  heat  flux  boundary  condition  obtained  is 
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In  addition  to  the  latent  heat  and  curvature  terms  on  the  left-hand  side  of  this  equation,  which 
were  also  present  in  equation  (72),  we  now  have  terms  involving  density  jumps  (e.g.  see  Burel- 
bach,  Bankoff  and  Davis  [36]).  The  last  term  on  the  right-hand  side  has  the  effect  of  reducing 
the  effective  latent  heat,  while  the  curvature  term  and  pressure  terms  may  contribute  positively 
or  negatively.  We  note  that  the  terms  on  the  right-hand  side  of  this  balance  also  appear  in  the 
modified  Clausius-Clapeyron  relation  (102)  and  so  can  be  related  to  the  viscous  dissipation  and 
the  interface  kinetics  effects. 
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5.3.4  Phase-Field  Equation 


The  leading-order  phase-field  equation  is  the  same  equation  (41)  as  in  the  equilibrium  analysis. 
An  analysis  of  the  phase-field  equation  may  be  pursued  in  a similar  manner  to  that  in  the 
equilibrium  case  in  Section  4.2  and  the  previous  nonequilibrium  case  in  Section  5.2.4.  An 
alternate  way  [1]  of  deriving  the  modified  Clausius-Clapeyron  relation,  however,  is  to  examine 
the  jump  in  Gibbs  free  energy  g.  One  expression  for  g follows  from  its  definition  in  (20)  and 
another  follows  directly  from  the  phase-field  equation  (31).  Equating  these  two  expressions  gives 
the  Clausius-Clapeyron  relation  for  a curved  interface  in  nonequilibrium,  as  we  show  below. 

Upon  expanding  the  definition  (20)  as  g — g^)  -I-  tg^  [parentheses  around  the  subscripts 
distinguish  between  the  terms  in  the  expansion  and  our  previously  defined  g0  in  equation  (7)] 
we  find,  using  equation  (85)  that  the  jump  in  p(0)  vanishes.  At  (9(e),  using  equations  (85),  (92) 
and  that  0o 1 1 = 0,  the  definition  of  g gives 
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where  we  have  used  r(0)  - 0 and  r(l)  — 1 but  otherwise  have  not  specified  the  form  of  r(4>). 

We  next  examine  the  jump  in  g as  obtained  directly  from  the  phase-field  equation  (31).  The 
leading-order  phase-field  equation  can  be  written 
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Since  g = g{0,p , 4>)  we  have 
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Noting  that  0O  is  independent  of  £ and  using  equations  (85)  and  (96)  we  obtain 
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This  gives  <?(o)ls  — 0 which  is  consistent  with  the  definition  of  g. 

At  the  next  order,  the  phase-field  equation  requires  that  g^  satisfy 
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We  obtain  from  equation  (97)  at  O(e)  that 
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where  we  have  used  that  0O  is  independent  of  £ and  that  dg(Q)/d6  = 0 [consistent  with  the 
dimensionless  entropy  in  equation  (19)  having  no  contribution  at  0(1)  when  (3  = 0 and  A = 
O(e)].  Some  additional  manipulations  of  (100)  show  that 
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Now  if  we  equate  the  two  results  (95)  and  (101)  we  find  the  modified  Clausius-Clapeyron 
relation  (in  dimensional  form)  that  includes  nonequilibrium  and  curvature  effects 
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These  integrals  are  evaluated  for  special  cases  in  Appendix  E. 

Notice  that  the  terms  in  the  first  line  of  equation  (102)  represent  effects  that  are  present  in  the 
equilibrium  case  while  those  in  the  second  line  represent  nonequilibrium  effects  associated  with 
the  change  of  phase.  Also  notice  that  M/lpi  = 0L/TMPmob  as  defined  in  (79).  As  a consequence 
of  our  quasi-incompressibility  assumption  [see  equation  (6)]  the  effects  of  compressibility  of  the 
bulk  fluid  are  not  included  in  this  equation.  We  refer  the  reader  to  Baldoni  and  Rajagopal 
[45],  where  the  effects  of  isothermal  compressibility  have  been  included  in  a modified  Clausius- 
Clapeyron  relation  for  a solid-liquid  system. 
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6.  Conclusions 


In  this  paper  we  have  examined  sharp-interface  limits  of  a phase-field  model  of  solidification  for  a 
pure  material  that  includes  the  effects  of  fluid  motion.  The  model,  given  by  equations  (29)-(32), 
was  developed  in  references  [1]  and  [22].  Here  we  considered  a simplified  version  that  retains 
the  square-gradient  and  double-well  terms  in  the  internal  energy  functional  but  excludes  them 
in  the  entropy  functional.  The  thermodynamically  consistent  derivation  given  in  [1]  and  the 
identification  of  the  distinguished  limits  which  lead  to  useful  free-boundary  problems  given  in 
this  paper  together  provide  a validated  phase-field  model  of  solidification  with  convection.  This 
paper  gives  the  first  sharp-interface  analysis  which  addresses  fluid  flow  and  solidification  in  a 
phase-field  model.  However,  we  note  that  a similar  sharp-interface  calculation  for  a binary  fluid 
(Cahn-Hilliard  model)  in  the  inviscid  case  has  been  carried  out  by  Lowengrub  and  Truskinovsky 
[16]- 

The  first  important  result  in  the  present  work  follows  from  the  sharp-interface  analysis  of 
the  equilibrium  version  of  the  model  in  which  we  considered  the  classical  limit  A = O(e)  as 
e — > 0.  Here  the  phase-field  model  leads  to  the  Young-Laplace  condition  (46a)  relating  the 
jump  in  pressure  between  the  bulk  phases  to  the  interface  curvature.  Also,  a Clausius-Clapeyron 
condition  (46b)  relates  the  temperature  of  the  system  to  the  interface  curvature  and,  when  the 
bulk  densities  are  different,  the  pressure  as  well.  Therefore,  even  when  there  is  no  flow,  a density 
difference  between  the  bulk  phases  can  affect  the  interfacial  temperature. 

In  the  first  of  two  distinguished  limits  for  the  nonequilibrium  case,  we  used  A = (A(e),  ps/ Pl~ 
1 = 0(e)  with  all  other  parameters  treated  as  0(1),  as  e 0.  The  hydrodynamic  boundary 
conditions  of  the  associated  free-boundary  problem  are  the  standard  ones  expected  for  a fluid- 
fluid  interface  undergoing  phase  transition  [namely,  the  normal  flux  boundary  condition  (57), 
the  condition  of  continuous  tangential  velocities  (64)  and  the  stress  balance  conditions  (69)]. 
It  is  important  to  point  out  that  the  condition  of  continuous  tangential  velocities  and  the 
stress  balances  both  follow  from  the  momentum  equation.  Viscosity  is  an  essential  element 
that  ensures  the  continuity  of  tangential  velocity  because,  in  the  sharp  interface  limit,  the 


28 


interface  cannot  support  unbounded  shear  stresses  and  strains  across  it.  Without  viscosity  we 
expect  that  the  sharp-interface  analysis  would  allow  a discontinuity  in  the  tangential  velocity 
across  the  interface.  The  thermal  interfacial  conditions  obtained  in  this  case  are  also  of  the 
familiar  types  [namely,  the  continuity  of  the  bulk  temperatures  at  the  interface,  an  expression 
for  the  interface  temperature  (78)  and  a condition  involving  the  jump  in  heat  flux  across  the 
interface  (72)].  However,  we  found  that  the  heat  flux  boundary  condition  involves  an  interface 
curvature  term  proportional  to  surface  energy  70.  This  term  is  often  not  included  in  standard 
free-boundary  formulations.  Its  origin  here  is  the  source  term  in  the  energy  equation  (32)  related 
to  the  gradient  energy  coefficient  e#  and  is  present  with  or  without  flow  and  density  differences. 
The  appearance  of  this  term  implies  physically  that  the  surface  energy  of  the  interface  affects 
the  balance  of  energy  at  the  interface.  Umantsev  and  Roitburd  [39]  found  that  the  more  general 
coefficient  of  the  curvature  term  is  y0  — Tdjo/dT.  We  surmise  that  had  we  retained  the  square 
gradient  and  double-well  terms  in  the  entropy  functional  we  would  have  recovered  this  more 
general  term  in  the  sharp-interface  limit.  Finally,  the  interface  temperature  in  this  case  is 
modified  from  its  reference  value  Tm  by  interface  curvature  and  attachment  kinetics. 

This  limit  allows  for  the  interpretation  of  the  solid  as  a highly  viscous  liquid  if  we  take  the 
subsequent  limit  Hs/hl  °°-  To  leading  order  in  hs/^l  the  corresponding  free-boundary 
problem  is  satisfied  by  rigid  body  motion  in  the  putative  solid  phase  while  the  velocity  in 
the  liquid  phase  satisfies  the  Navier-Stokes  equations  subject  to  normal  mass  flux  and  no-slip 
conditions  at  the  interface  (see  the  discussion  at  the  end  of  Section  5.2.2). 

A second  distinguished  limit,  given  by  A = O(e)  and  Pr/ 7 = O(c),  is  also  possible  in  the 
nonequilibrium  setting.  This  limiting  case  may  be  more  appropriate  for  a fluid-fluid  or  fluid- 
vapor  system  where  viscous  effects  may  be  of  secondary  importance.  Here  the  hydrodynamic 
conditions  are  similar  to  the  previous  nonequilibrium  case  with  the  exception  that  the  viscous 
terms  are  not  present  in  the  stress  jump  condition  at  leading  order  [see  equations  (83),  (91) 
and  (93)].  The  condition  of  continuous  tangential  velocities  is  still  obtained,  reflecting  the  limit 
of  small  but  nonvanishing  viscosity.  The  thermal  interfacial  conditions  have  the  same  general 
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structure  as  those  in  the  previous  nonequilibrium  case  [see  equations  (94)  and  (102)]  but  now 
include  additional  effects  that  can  be  attributed  to  density  differences  between  the  solid  and 
liquid  phases. 

The  phase-field  model  developed  in  [1]  is  a so-called  generalised  phase- field  model,  as  it 
includes  additional  square-gradient  and  double-well  contributions  to  both  the  internal  energy 
and  entropy  functionals  of  the  system.  Similar  models  have  been  developed  by  other  authors 
[14,  33,  40,  46,  47].  McFadden,  Wheeler  and  Anderson  [30]  considered  various  sharp-interface 
limits  of  a generalised  phase- field  model  and  investigated  the  so-called  thin  interface  limit,  in 
which  A = 0(  1)  and  X0  = O(e)  as  e — » 0.  They  found  that  the  generalised  phase- field  models 
can  give  rise  to  additional  nonstandard  terms  in  the  thermal  interfacial  boundary  conditions. 
Here  we  examined  different  asymptotic  limits  both  with  A — O(e)  as  e -»  0 and  also  found 
that  the  leading  order  thermal  interfacial  conditions  involve  new,  albeit  different,  terms  directly 
related  to  the  square  gradient  and  double-well  terms  in  the  energy  functional. 
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Appendix  A.  Notes  on  Differential  Geometry 

We  employ  body  fitted  coordinates  (CSC^C3)  t°  the  interface  [4>  = 1/2  given  by  r = ro(C\C2,£)] 
(e.g.  see  Aris  [48]).  Here  C1^2  measure  distance  along  the  principal  directions  of  the  interface 
and  C3  measures  distance  in  the  direction  normal  to  the  interface.  Hence 

r = r0(CX2,t)  + C3n, 
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where  n is  the  unit  normal  to  the  the  interface.  The  corresponding  basis  vectors  are 

e*i  = (1  — C3^i)^i5  e*2  — (1  ~ C3^)^  e*3  — fb 

where  ti,t2  are  the  unit  tangent  vectors  to  the  interface  parallel  to  the  principal  directions. 
The  corresponding  metric  tensor  has  nonzero  components  g11  = (1  — C3 ^i)-2 , g22  = (1  — 
C3JC2)-2,g33  = l. 

The  gradient  of  </>  may  be  expressed  as 

d(f) 
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where  e k • e}  = 6j.  Hence  we  find  that 
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The  contravariant  components  of  V0  <g>  V</>  are 
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We  note  that  the  contravariant  components  of  the  unit  tensor  in  these  coordinates  is  given  by 

I ~ diag(gn , g22 , g33).  Hence  the  expansion  of  cr^  = ye  [||V</>|2/  — V0C8)  V</>]  in  the  inner 
region  (with  £3  = e£)  is  <x^  = e-1  <Tq  + ecrf  + (9(e2)j  where 
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In  a similar  way  the  expansion  of  crp  — - (p  + p*)  I in  the  inner  region  is  given  by  crp  — 


--i  \~p 


+ tcr\  + (9(e2)],  where  the  contravariant  components  of  crfi  and  cr\  are 
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and 
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The  viscous  part  of  the  stress  tensor  in  the  inner  region  is  likewise  given  by  r = e_1  [r0  + erx  + 0(e2)\, 


where  the  contravariant  components  of  r0  are 
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The  only  component  of  rx  needed  for  the  present  calculations  is  t33  which  is  given  by 
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We  note  that  the  components  of  the  divergence  of  a second  rank  tensor  A are  given  by 
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where  g = gn922933  = (1  — £3^i)2(1  — £3Ai2)2  and  are  the  so-called  Christoffel  symbols.  We 
now  expand  V • A in  the  inner  region.  We  note  that  the  only  Christoffel  symbols  with  a non-zero 
order  one  contribution  are 
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Hence  if  A = A0  + cAx  + 0(e2)  then 
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Appendix  B.  Surface  Energy 


The  surface  energy  7 may  be  defined  as  the  excess  Kramers  potential  (defined  by  f — g,  where  / 
and  g are  the  Helmholtz  and  Gibbs  free  energies  per  unit  mass)  of  an  interface  at  the  equilibrium 
temperature  Tm  at  pressure  pF.  If  we  introduce  a dividing  surface  at  2 = zF  the  surface  energy 
is  given  by 
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which  is  independent  of  the  position  of  the  dividing  surface.  In  equilibrium  the  the  stress  in  the 
interface  is  zero  and  hence  m — 0 and  thus  from  (36) 
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and  hence  we  may  express  the  surface  energy  as 

r+o o r+00  /7<bn\  2 

7 = y_oo  e2F{TM)(t>]dz  = 670  dQ,  (B3) 

where  70  = \€-f{Tm)\J Pl/ (2a)  corresponds  to  the  surface  energy  when  p$  = Pl  [see  equa- 
tion (11)]. 


Appendix  C.  Energy  Equation  Analysis:  Ps  I PL- 1 = 0{e) 

Here  we  outline  the  calculation  for  the  derivation  in  the  sharp-interface  limit  of  the  heat  flux 
boundary  condition  for  the  case  A = Ac  and  Ps/pl  ~ 1 = O(e)  with  all  other  parameters  0(1) 
as  e — » 0.  In  this  case,  we  note  that  v = S = S/X  = ( S/A)e~ 1 and  from  equation  (18)  that 
e = 
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At  0(e~2)  the  energy  equation  (32)  gives  that  the  leading  order  temperature  in  the  interfacial 
region  is  independent  of  £ and  hence  may  be  expressed  as  0o(CSC2A)- 
At  the  next  order  we  find  that  the  energy  equation  gives  that 
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and  that  the  continuity  equation  gives 
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Equation  (Cl)  may  be  simplified  by  using  (C2),  the  leading  order  phase-field  equation  (74),  as 
well  as  the  expression  for  P0  given  by  (63)  to  yield 
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We  integrate  this  once  and  apply  the  matching  conditions 
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Noting  that  the  integral  in  this  expression  evaluates  to  1/6,  the  dimensional  form  of  this  interface 
condition  is 


kLh  • V7l  — ksh  ■ VTs  — (w0|s  * A — Vn)  [plL  + 7o(/Ci  + /C2)j . 


(C6) 


34 


Appendix  D.  Energy  Equation  Analysis:  Ca  = 0(e) 

Here  we  outline  the  calculation  for  the  derivation  in  the  sharp-interface  limit  of  the  thermal 
boundary  conditions  for  the  case  Pr/j  = Ce,  X — Ae  with  all  other  parameters  0(1)  as  e — > 0. 
As  before,  we  note  that  u = S = S/X  — ( S/A)e~l  and  from  equation  (18)  that  e — C9(e-1). 

As  for  the  previous  case  with  nearly  matched  densities,  at  0(e~2)  the  energy  equation  (32) 
gives  that  the  leading  order  temperature  in  the  interfacial  region  is  independent  of  £ and  hence 
may  be  expressed  as  0o(C1>C2^)- 

At  next  order  C9(e_1)  the  energy  equation  (32)  gives  that 
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where  Dx  represents  the  G(e  *)  contribution  of  the  dissipation  terms  e|-^V2<f)  + 'H  and  is  given 
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Here  Px  is  given  by  equation  (92)  and  C3(£1,  (2,  t)  = pols  + p*  + Jo  u0|s  • n.  The  normal  velocity 
correction  Uf.  which  can  be  obtained  by  examining  the  continuity  equation,  is  given  by 
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Using  this  result  we  can  integrate  equation  (Dl)  once  and  apply  the  matching  conditions  (C4) 
to  obtain 
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Next,  we  use  the  result  for  Uf  to  simplify  Di  and  find 
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If  we  note  that  C/q  - Uq\s  = Jo/p($o)  — PlJq/ps  it  follows  that 
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If  we  integrate  the  last  term  in  equation  (D5)  by  parts,  use  equations  (D4)  and  (D7)  and  again 
note  that  (a<E>0/a£)2  = p(^o)Hrn(^0)  we  find  that 
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In  dimensional  form  this  heat  balance  is 
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Here  we  have  assumed  that  p($0)Hm($o)  = <f>o(l  — d>o)2  so  that  the  integral  in  equation  (D8) 
evaluates  to  1/6. 
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Appendix  E.  Integrals  in  Clausius-Clapeyron  Relation 


Suppose  p{<f>)Hm{<f>)  = 02(  1 - 4 >)2  and  p(0)  = 1 + {ps/pL  ~ 1 M0).  Then  if  r{<j>)  = 02( 3 - 2 <f>) 

we  have 
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and  if  r(4>)  = 4>  with  Si  = ( Ps/pL  - 1)  we  have 
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If  r(4>)  = 4>  and  p{4>)  = 1 + S2r{4>)  where  S2  — ( Ps/pl  ~ 1)  then 
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If  we  fix  S2,  Z2  increases  as  Si  increases  or  decreases  away  from  zero.  We  note  that  for  the 
expression  in  equation  (E3),  ^ ~ \s\  (1  + s2/2)  when  sj  < 1.  A similar  behavior  is  found  for 
Z2  with  r(4>)  — 02( 3 - 20). 
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